Finite Element Analysis of Trabecular Damage Adjacent to Interbody Fusion Devices

INTRODUCTION: Surgical intervention for the treatment of low back pain over the last decade has often included the use of interbody fusion devices (IFDs). Implanted into the intervertebral space, these devices serve as structural support until fusion occurs. Numerous experimental investigations have evaluated the acute kinematic stability and/or monotonic strength of the implant-vertebrae constructs. A few finite element analyses (FEA) have focused upon bone-implant interactions; however none have incorporated complex material behaviors even though this interface directly affects postoperative stability and therefore success. The objectives of this investigation were to determine the difference in predicted damage accumulation adjacent to two different IFD geometries and to determine the effect of the inclusion of viscous material behavior on the model predictions relative to results from a previous experimental and histological study [7]. METHODS: Generic cylindrical and hexahedral IFD geometries were incorporated into idealized 2D, plane strain models generated from QCT data of the experimental specimens (Fig 1). Medial-lateral and inferior-superior symmetry were utilized for computational efficiency. The vertebral endplates and were respectively assumed to have uniform thicknesses of 0.3 and 0.5mm. The four materials identified within each model (vertebral shell, vertebral centrum, vertebral endplate and implant) were assumed to be homogeneous and isotropic. Elastic modulus and Poisson’s ratio values from the literature were assigned to the vertebral endplates, the vertebral shell and the cylindrical implant (1000MPa, 0.2; 5000MPa, 0.3, and 110GPa, 0.3 respectively). The central cavity of the hexahedral implant was modeled as a solid cross-section with a reduced modulus (2.9GPa) to compensate for this structural modification. The elastic properties for the cancellous vertebral centrum were based upon QCT measurements and a published modulus-density relationship [4] and ranged from 240365MPa; a Poisson’s ratio of 0.2 was assigned. The inelastic cancellous material behaviors included plasticity with linear hardening (PLH) and two-layer viscoplasticity (2VP) (ABAQUS 6.5). The isotropic yield stress (σy) was determined using relationships found in the literature [3] and the hardening modulus (H) determined as the slope of the line between the ultimate and yield stress-strain values for each specimen. In addition to the plasticity parameters, the 2VP material behavior (Fig 2) requires a creep law for the viscosity coefficient (c), the elastic moduli (Ev, Ep and Ei [=Ev+Ep]) and time constant (τ). τ was determined from the fit of a five parameter linear viscoelastic model to the experimental data and (1-Ev/Ei) chosen as 0.25 (2VP) based upon prior work in our laboratory. Using a time-hardening power law, absent the transient and non-linear viscoelastic terms, the viscosity coefficient was then determined for each model. The bone-implant interface was modeled as a contact pair with finite sliding. The tangential behavior was implemented with an isotropic friction coefficient of 0.2 [5; 6]. A load-hold-unload compressive displacement cycle was applied to each model via the endplate nodes. Hold displacement amplitudes were chosen to generate either 1.0% or 2.5% average strain. Eight analyses were executed; two cancellous material behaviors (PLH & 2VP) for the two IFD geometries to two displacements. Structural force and displacement values were obtained from the endplate nodes. Because an explicit material model for damage accumulation in cancellous bone was not available, it was assumed that inelastic strain represented concomitant plastic flow and damage accumulation. Accordingly, the maximum plastic equivalent strain (PEEQ) for each model was obtained from its corresponding integration point. RESULTS: For both devices, the predicted damage zone was localized near the bone-implant interface. In the PLH material model at both 1% and 2.5% structural strain, the typical PEEQ strain distribution in the cylindrical IFD model fully encompassed the cancellous region below the device. The distribution for the hexahedral device showed little PEEQ strain immediately below the device and a maximum at the lateral corner (Fig 3). On the other hand, the force-displacement curves showed little difference between the two IFD geometries at either 1 or 2.5%.